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INTRODUCTION 


A.  BACKGROUND 

1.  The  Importance  of  Stealth 

The  concept  of  stealth  has  been  around  since  long  before  man  made  his  appearance 
on  earth.  Nature  realized  that  the  survival  of  its  creatures  depended  on  their  blending  in 
with  the  background.  Evolution  has  provided  the  means  for  nature’s  creatures  to  develop 
and  refine  stealth  capabilities.  It  has  been  during  the  last  fifty  or  so  years  that  the  military 
forces  of  the  world  have  taken  a  keen  interest  in  stealth  [1],  The  results  were  displayed 
by  the  success  of  the  FI  17A  stealth  aircraft  during  the  Persian  Gulf  War  [2,  3,  4]. 

The  success  of  these  aircraft  during  the  Gulf  War  can  be  attributed  to  many 
factors.  By  far  the  most  important  factor  was  possession  and  use  of  stealth  technology. 
The  U.S.  military  had  it  and  the  enemy  did  not.  More  importantly,  the  enemy  did  not 
possess  the  technology  to  defeat  a  stealthy  airborne  platform.  That  was  in  1991,  more 
than  six  years  ago.  Since  then  stealth  technology  has  been  ever  so  slowly  let  out  of  its 
“black  box.”  It  is  no  longer  “our  own  little  secret  weapon”  that  only  the  U.S.  military 
possesses,  and  there  is  no  guarantee  that  during  the  next  conflict  our  adversaries  will  not 
possess  anti-stealth  technology.  In  fact  ultra-wideband  radar  (UWB)  was  being 
developed  for  detection  of  stealth  platforms  before  the  Gulf  War  began  [5],  In  1994  the 
United  States  Air  Force  admitted  that  some  radars,  including  some  mobile  units,  could 
detect  the  B-2  bomber  [6],  And  as  one  would  expect,  the  Russians  are  also  developing 
their  own  anti-stealth  technology  [6],  The  bottom  line  is  that  if  stealth  platforms  are  to 
remain  stealthy,  they  need  to  continue  to  decrease  their  observability. 
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The  idea  of  stealth  is  comprised  of  several  aspects  of  observability  including 
infrared  radiation,  optical,  acoustic,  and  radar  echo.  The  objective  is  to  minimize  the 
observables  associated  with  your  platform  so  as  to  reduce  the  chances  of  being  detected 
by  the  enemy.  While  current  stealth  platforms  have  incorporated  reductions  in  each  of 
these  areas,  the  most  emphasis  has  been  placed  on  reduction  of  radar  signature.  This  is  a 
consequence  of  radar  being  the  primary  means  of  detection  for  the  U.S.  military,  its  allies, 
and  its  adversaries.  That  being  the  case,  this  thesis  attempts  to  take  one  small  step  in  the 
direction  of  radar  cross  section  reduction  by  determining  the  induced  source  distribution 
on  a  body  which  produces  scattered  electromagnetic  radiation. 

2.  Goals  of  the  Research 

Radar  cross-section  (RCS)  reduction  is  a  key  element  of  low-observability. 

Current  techniques  used  in  determining  the  RCS  of  a  platform  rely  on  analysis  of 
measurements  made  in  the  far-field.  Generally  these  measurements  provide  a  gross  picture 
of  the  platform’s  overall  RCS  as  a  function  of  viewing  angle.  This  information  enables 
engineers  to  then  modify  the  design  in  an  attempt  to  reduce  its  RCS.  This  becomes  an 
iterative  process  with  design  modifications  leading  to  measurements  which  lead  to  more 
design  modifications. 

At  the  other  end  of  the  spectrum,  measurements  made  using  a  probe  on  or  near  the 
outer  surface  of  the  body  may  influence  the  quantities  being  measured.  In  other  words  the 
act  of  measuring  would  induce  errors  into  the  quantities  being  sought. 

The  focus  of  this  study  is  to  more  accurately  determine  the  source  of  scattering  on 
a  body.  This  will  be  done  by  evaluating  the  viability  of  back-propagating  measurements 
to  an  axisymmetric  body  to  image  the  source  of  scattering  on  the  body.  This  analysis  will 
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enable  more  precise  location  of  scattering  sources  than  does  gross  analysis  of  the  far-field. 
Once  the  local  scatterers  have  been  identified  it  is  then  theoretically  possible  to  remove  or 
reduce  them  in  the  overall  effort  of  radar  cross-section  reduction. 

The  objective  of  this  study  is  to  investigate  the  viability  of  back-propagating 
electromagnetic  field  measurements  to  an  axisymmetric  body  in  order  to  determine  the 
source  distribution  for  radiated  power  from  the  body.  If  this  is  indeed  possible,  then  the 
follow-on  objectives  include  determining  the  range  and  spatial  resolution  at  which 
measurements  must  be  made  in  order  to  provide  meaningful  results. 


B.  SOURCE  IMAGING 

1.  Acoustics  Work  with  Cylinders 

A  brief  history  of  this  subject  starts  with  researchers  attempting  to  determine 
sources  of  acoustic  noise  on  a  body.  It  was  shown  that  acoustic  intensity  could  be  used  to 
localize  sources  of  sound  on  structures  which  radiate  to  the  far  field.  A  new  quantity 
called  the  supersonic  acoustic  intensity  vector  was  defined  and  its  application  to 
measurements  on  plate  and  cylinder-like  structures  was  demonstrated  [7], 

“Supersonic  intensity  is  composed  only  of  wave  components  which 
radiate  to  the  far  field  (supersonic),  with  the  non-radiating  (subsonic) 
components  eliminated.  The  normal  component  of  this  supersonic  intensity 
vector,  measured  in  the  extreme  near  field  or  on  the  surface  of  the 
structure,  provides  an  accurate  tool  for  locating  regions  (“hot  spots”)  on 
the  structure  which  radiate  to  the  far  field.  Furthermore,  the  supersonic 
intensity  provides  an  accurate  quantification  of  these  source  regions, 
providing  a  ranking  of  the  strength  of  the  identified  source  regions  as  a 
function  of  frequency.  This  identification  and  ranking  provides  a  powerful 
new  tool  in  the  understanding  and  control  of  radiated  noise.”  [7] 


Coupled  with  the  finding  that  acoustic  source  regions  could  be  determined  through 
back-propagation  was  that  the  results  were  obtained  in  two  dimensions.  Using  a  cylinder 
as  the  scattering  body,  scattering  data  measurements  were  made  on  an  arc  of  a  circle  with 
diameter  on  the  axis  of  the  cylinder  [8].  In  the  spherical  ( r,0,<f> )  coordinate  system  this 
translates  to  measurements  over  6  at  a  constant  <}> .  The  restricted  case  of  axisymmetric 
fields  and  currents  is  the  next  level  of  difficulty  before  working  with  full  three  dimensional 
shapes. 

2.  Extension  to  Electromagnetics  for  Cylinder  Geometry  Case 

The  discoveries  employing  supersonic  modes  to  enhance  acoustic  imaging  on 
cylindrical  shells  has  been  translated  into  the  realm  of  electromagnetics  for  the  case  of 
cylindrical  geometry.  It  was  shown  that  superluminal  modes  which  satisfy  |k2|  <  k  and 
have  faster-than-light  axial  propagation  velocity  provide  time-average  power  to  the  far 
field.  An  optimal  propagation  deconvolution  filter  was  implemented  to  remove  the  effects 
of  the  subluminal  modes  during  back  propagation.  [9] 

3.  Extension  to  General  Non-separable  Geometries  Using  Finite  Element 

Methods 

The  goal  of  this  research  is  to  extend  the  concepts  developed  in  previously 
conducted  research  to  bodies  with  non-separable  geometries.  By  definition  these  bodies 
have  shapes  that  do  not  lend  themselves  to  explicit  rigorous  solutions  using  separable 
modes,  such  as  cylindrical  or  spherical  wave  functions.  A  finite  element  method  (FEM)  is 
used  to  relate  the  conditions  at  the  boundary  (measuring  region)  to  conditions  on  the 
surface  of  the  body,  in  this  case  relating  the  scattered  field  to  the  source  current  that 


generated  it.  The  number  of  elements  used  is  chosen  so  that  the  contour  of  the  body  is 
closely  matched.  Figure  1  shows  a  meridian  plane  of  a  sphere  surrounded  by  a  finite 
element  mesh  which  connects  the  surface  of  the  sphere  to  the  arc  of  a  circle  with  diameter 
on  the  axis  of  the  sphere,  as  previously  mentioned.  In  Figure  2,  the  body  represented  is  a 
cone.  The  more  contours  on  the  body,  the  greater  the  number  of  elements  required  to 
closely  match  those  contours.  Ultimately  the  accuracy  of  the  solution  is  determined  by 
the  number  of  elements  used. 
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Figure  1 .  Sphere  (Meridian  Plane)  Surrounded  by  Finite  Element  Mesh 
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Figure  2.  Cone  (Meridian  Plane)  Surrounded  by  Finite  Element  Mesh 


In  the  chapters  that  follow,  the  source  current  on  the  surface  of  an  axisymmetric 
body  is  determined  through  back-propagation  of  the  scattered  field  it  generates.  In 
Chapter  II,  the  scattered  field  (measured  field)  is  determined  through  numerical  integration 
of  a  surface  current  induced  on  an  axisymmetric  body  by  a  locally  placed  elemental  dipole. 
Chapter  III  details  the  back-propagation  of  the  scattered  field.  Chapter  IV  provides  an 
analysis  of  both  the  numerical  integration  results  and  the  back-propagation  results. 
Conclusions  are  presented  in  Chapter  V.  An  appendix  which  includes  all  pertinent 
computer  code  is  also  included  to  assist  in  the  understanding  and  continuation  of  this 
research. 
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H.  GENERATION  OF  SCATTERED  FIELD  (MEASURED  FIELD) 


A.  BACKGROUND 

Before  the  scattered  field  can  be  back-propagated  it  must  be  measured.  To 
simulate  the  process  and  control  the  sources  of  error,  the  measured  field  is  computed  to  a 
specified  level  of  accuracy.  While  reducing  potential  error  associated  with  field 
measurements,  generating  fields  leads  to  an  additional  hurdle,  namely  how  to  generate  a 
scattered  field  from  an  arbitrary  axisymmetric  body.  This  problem  is  solved  by  inducing  a 
surface  current  on  the  body  and  integrating  it  to  find  the  resulting  scattered  field. 

Potential  error  is  introduced  as  a  result  of  integration  accuracy.  A  sphere  is  first  chosen  as 
the  arbitrary  axisymmetric  body  in  order  to  test  the  accuracy  of  the  general  source 
integration  algorithm. 

A  particular  scattered  field  is  generated  by  placing  a  radial  directed  elemental 
dipole  in  the  vicinity  of  the  metal  sphere.  The  dipole  induces  a  surface  current  Js  on  the 
sphere  which  in  turn  generates  scattered  fields  Es  and  Hs .  The  measured  field  H  at 
some  specified  distance  from  the  sphere  is  composed  of  the  field  scattered  (Hs )  from  the 
sphere,  and  the  field  incident  (Hj )  from  the  elemental  dipole.  The  ‘exact’  induced 
currents  and  scattered  fields  can  be  computed  for  this  test  case  and  used  to  validate  the 

accuracy  of  the  numerical  integration  and  finite  element  algorithms.  Calculation  of  the  H 
field  is  considered  in  Section  C.  It  is  generated  using  the  program  DIPSPHR2.M  found  in 
the  Appendix.  Field  integrations  for  axisymmetric  currents  are  detailed  in  Section  D  and 
integration  results  are  provided  in  Section  E  of  this  chapter. 
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B.  SPECIALIZATION  TO  AXISYMMETRIC  FIELD  CASES 

Determining  the  source  distribution  on  bodies  with  non-separable  geometries  is  the 
ultimate  objective  of  research  in  this  area.  This  thesis  reduces  the  problem  from  three 
dimensions  to  two  dimensions  by  assuming  axisymmetric  fields  generated  by  axisymmetric 

currents  on  a  body  of  revolution.  There  are  two  special  cases  to  be  considered:  the  TE , 

— — 

case  in  which  the  E  field  and  the  surface  current  are  transverse  to  <f> ,  and  its  dual,  the 

-  A  A 

TM+  case  in  which  the  H  field  is  transverse  to  <f>  and  the  surface  current  is  <j>  -directed. 
These  cases  are  depicted  in  Figure  3. 


Restricted  Case:  TEt  Dual:  TM+ 

Figure  3.  Axisymmetric  Field  Cases 


This  thesis  will  investigate  the  TE+  case.  The  fields  for  this  special  case  are. 


H(j>,z)  =  H+(p,z)<f> 

(la) 

E(p,z)  =  Ep(p,z)p+E,(p,z)z 

(lb) 
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C.  OFFSET  DIPOLE  TO  GENERATE  SURFACE  CURRENT  ON  SPHERE 
A  radial  directed  dipole  located  in  the  vicinity  of  a  metallic  sphere  is  used  to 
generate  surface  currents  on  the  sphere.  These  surface  currents  are  then  integrated  to  find 
the  scattered  field  at  a  specified  distance.  The  dipole  case  is  also  used  to  test  the 
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integration  accuracy  by  providing  a  nearly  exact  scattered  field.  Figure  4  shows  an 
elemental  dipole  located  in  the  vicinity  of  a  metal  sphere  in  the  spherical  coordinate 


system. 


Figure  4.  Dipole  in  Vicinity  of  Metal  Sphere 

The  procedure  for  determining  the  surface  current  consists  of  several  steps.  The 
Ei  and//,  fields  due  to  the  dipole  are  represented  in  ( r,0,(j> )  coordinates  of  the  centered 
sphere.  The  scattered  fields  Es  and  Hs  due  to  the  induced  currents  on  the  sphere  are 
then  expanded  using  spherical  harmonic  expansions  with  unknown  coefficients.  Enforcing 


10 


rx(Ei+Es)  =  0  =>  Eton  =  0  on  the  sphere  surface  allows  solution  for  the  expansion 
coefficients.  Finally  Js=rx  (Hi  +  Hs)  is  used  to  find  induced  surface  currents  which 

produce  Es  and  Hs  fields  scattered  from  the  sphere. 

An  elemental  dipole  placed  at  the  origin  produces  the  following  fields  [10]: 


E 


rd 


Vpldl_ 

In 


cos  0d 


(6a) 


riJdl 


Ee  = 

“  4  n 


r>  Jr 

sin  9,  - +  —5-  - 


H 


Idl_ 
4  n 


sin  9d 


(6c) 


The  position  of  the  dipole  must  then  be  translated  to  the  position  z  =  z0 .  This 


translation  is  accomplished  using  the  following  transformation  formulas: 


rd-  yjr2  +z20  -  2 rz0  cos  0 

(7a) 

„  rcos0-zo 
cos9d  - 

r4 

(7b) 

rsin# 

sin#,  - 

rd 

(7c) 

EB  =  Egd  cos(0 -0d)~  Erd  sin(<9 -  0d) 

(8a) 

Er  =  Ee<j  sin(0 -9d)  +  Er<t  cos {6-  0d) 

(8b) 

(8c) 
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The  transformation  formulas  in  Equations  (7)  and  (8)  are  then  used  with 
Equation  (6)  to  express  the  fields  due  to  the  dipole  on  the  surface  of  the  metal  sphere, 

Eei(a,0 ) ,  Eri(a,&) and  .  This  field  will  induce  a  surface  current  J  =  Jg{6)§ 

on  the  surface  of  the  sphere  which  will,  in  turn,  generate  scattered  fields  Ees,  Ers,  and 
H+s.  The  tangential  E  field  at  the  sphere  surface  will  exactly  cancel  the  tangential  E 

field  due  to  the  dipole  thus  producing  a  total  tangential  E  field  of  zero 

E„,(a,e)  =  -E,l(a,6)  (9a) 

The  surface  current  will  also  satisfy  J  =  nx  Huui  on  the  surface: 

j  =  ?x  =  “  <**  +  *  ft  = 

=5  J,(0)  =  (9b) 

This  Jg(0)  is  what  needs  to  be  integrated  to  find  the  radiated  H^s(r,8)  at  locations  off 
the  surface  of  the  metal  sphere.  The  result  is  H^s  which  will  be  used  to  compare  results  of 

the  integration  covered  in  the  following  section. 

Note  that  although  Je  is  produced  by  the  total  on  the  metal  surface,  it 

generates  only  the  scattered  field.  The  scattered  field  has  an  Ees  which  exactly  cancels 
out  the  Eei  of  the  incident  dipole  field  on  the  sphere  surface. 

To  determine  the  “exact”  solution,  the  scattered  fields  can  be  expressed  as 
weighted  sums  of  spherical  harmonics  [11], 

t 

*  H{2)  (Or) 

E„(r,0)  =  MXa.  ’  (10a) 

n=l  P' 
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H„(r,e) = ±a.22jP-r>.<) cos#)  (10b) 

n= 1  P* 

where  H^ifir)  =  -  jYn  (ftr)  is  the  “Riccati”  spherical  Hankel  function  of  the 

second  kind  representing  outbound  waves.  Pj(cos0 )  is  the  “associated  Legendre 
function”  and  is  the  m  =  1  case  of  P”( cos0) .  To  numerically  solve  the  problem, 
truncate  the  series  in  Equation  (10)  and  solve  for  N  values  of  complex  an .  Using 


Equation  (9a)  and  defining  Dn  = 


pa 


gives 


jThT.anDnPl(?os0)  =  -Eei(a,0)  (1 1) 

n= 1 

The  an  are  found  by  use  of  orthogonality  of  the  associated  Legendre  functions, 


P„m(cos0) .  For  the  special  case  m  =  1 , 


|  P„’  (cos  0)P}  (cos0)  sin  0d0 
0 


f  2n(n  + 1) 


i  2n  +  l 
I  0  l*n 


=  n 


We  can  sift  out  coefficients: 

=  -]E„Xa,ff)P:(ooSe)mm 

2n+l  J 


(  h  (2w  + 1) 

*  J'hD.  2n(n  +  l) 


(12) 


(13a) 


(13b) 


The  /n  integration  can  be  performed  numerically  to  any  accuracy  since  E6j(a,&)  and 


/('(cos#)  can  be  obtained  using  as  many  # -points  as  desired.  Once  the  an  are  found  the 
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Jg(6)  can  be  computed  at  any  0  value  by  use  of  Equation  (9b)  with  H+s  given  by 
Equation  (10b). 

D.  FIELD  INTEGRATION  FOR  AXISYMMETRIC  CURRENTS 

As  noted  in  the  introduction,  the  first  step  is  to  measure  the  field  that  will  be  back- 
propagated  to  the  surface  of  the  axisymmetric  body.  For  our  purposes  we  will  assume  an 
axisymmetric  surface  current  distribution  on  the  body.  The  “measured”  fields  will  actually 
be  determined  through  the  integration  of  the  surface  current.  In  this  context  the  scattered 
field  is  generated  through  integration  rather  than  through  illumination  and  measurement. 

The  surface  of  an  arbitrary  axisymmetric  body  of  revolution  can  be  defined  by  use 
of  the  generating  contour  p’(z')  using  the  ip,<j>,z)  cylindrical  coordinate  system  as 
shown  in  Figure  5.  Source  points  on  the  body  are  defined  by  individual  (p,z)  pairs.  The 
field  point  in  space  is  located  at  (p,<p,z) .  The  axisymmetric  surface  current  distribution  is 
Js(r')  =  Jt{p',z')i  +  J ,z')<f> .  The  unit  vectors  f  and  <j>  are  tangent  to  the  surface  at 
each  point  and  both  Jt  and  are  constant  with  changing  <t> . 
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The  fields  will  be  numerically  evaluated  via  the  vector  potential  formulation  given 
by  [11]: 


H(r)  =  —  V  x  A(r ) 


(14a) 


E 


— VxH  =  jcoA- 


V(V-A) 


1 


(14b) 


with  R-r-r'  and  R  -  \r - r'\ .  For  field  points  and  source  points  with  coordinates 
( p,<t>,z )  and  ')  the  law  of  cosines  gives 

R  =  -Jp2  +  p'1  -  ~2pp'  cos(^  -  ^')  +  (z  -  z')2  (16) 

The  curl  operator  in  Equation  (14a)  is  taken  inside  the  integral  in  Equation  (15)  and  note 

-jpR  /  i—  - » 

that  it  only  operates  on  the  e  term  (unprimed  position  in  R  =  \r  -  r'\  ) 

H<r)  =  J/  vl^-J  x  J,(?)dS'  (17a) 

surface  ^  ' 


where 


-R 

fl+  JfiR) 

R  J 

l  R3  J 

-/A* 


(17b) 


Thus  the  general  solution  is: 

1 


=  JJ[j,(r>)x  R][-  +jfR  1  e'i/utdSl 


An 


surface 


R 


(18) 


The  next  step  is  to  derive  the  explicit  integrations  to  be  performed.  To  do  so  we 
select  r  in  the  x-z  plane,  with  </>  =  0  as  shown  in  Figure  6.  Since  Js  and  the  resulting 
fields  so  derived  are  axisymmetric  (not  functions  of^)  we  lose  no  generality  by  doing  this. 
In  fact,  if  we  find  Hx,Hy  and  Hz  at  this  field  point  we  note  that  for  other  locations  where 
<f)  *  0  the  substitutions  Hx  — »  Hp,Hy  — »  H+  and  Hz  — >  Hz  apply.  Also  note  that 


r  =  px  +  zz 

r'  =  p'cos^'x  +  p'  sin  <f>'y  +  z’z 


(19a) 

(19b) 
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R  =  r-r'  =  (p- p'cos0')x  -  p'  smfiy  +  (z  -  z')z 


(19c) 


♦ 


Figure  6.  Geometry  for  Explicit  Integrations 
The  field  integrations  will  be  performed  numerically  by  breaking  the  axisymmetric 
surface  into  a  large  number  of  flat  circular  rings  as  shown  in  Figure  7a.  These  rings  are 
stacked  to  approximate  the  surface  of  revolution.  Each  ring  is  “flat”  in  the  direction  of 

i  and  curved  in  the  direction  of  ^  as  depicted  in  Figure  7b. 
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Figure  7a:  Stacked  Circular  Rings  form  Axisymmetric  Body 


z 


Figure  7.  Decomposition  of  Axisymmetric  Body 
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Consider  now  evaluating  H(r)  in  Equation  (18)  due  only  to  Js(r')  =  Jt  (s)t  over 
a  single  ring  as  in  Figure  7b.  Over  the  ring  we  will  approximate  J,(s)  as  piecewise  linear, 
specified  using  its  values  at  the  top  and  bottom  of  the  ring,  Jt  (5, )  and  Jt  (s2 )  via 


(s-s.)  (s-s.) 

J,  (s)  =  J,  ( s,)f - ^  +  J,(s2  )j—~\ 

v5l  S2j  \S2  S\) 


(20) 


This  equation  is  a  straight  line  interpolating  function  between  known  values  at  5,  and  s2 . 

A  single  ring  is  shown  in  Figure  8a  where  s  is  the  path  length  over  the  ringed  surface  from 
the  +  z  axis.  The  piecewise  linear  variation  of  J,  ( s  )  between  s,  and  s2  is  shown  in 

Figure  8b. 

To  evaluate  Equation  (18)  for  Js(r')  =  J,(s)t  over  a  single  ring,  we  note  that 


t  =  cosacos0'x  +  cosasin^'j)-sina£ 


(21) 


where  cos  a  =  t  •  p  and  -  sin  a  =  i  •  z .  After  simplifying, 

txR  =  x(tyR t  -  t2Ry)+ y(tzRx  -  txRz)  +  z[txRy  -  tyRx)  (22) 

and  substituting  Equation  (22)  into  Equation  (18)  with  restriction  to  one  ring  gives, 


1 

Hx  =  —  J  Jt  (-s)[(z  -  z ')  cos  a  -  p'  sin  a\Fs{s)p'ds 


(23a) 


t  s2 

Hy  =  —  J  J,  (5){[p'  sin  a  -  (r  -  z')  cos  a\Fc  (5')  -  p  sin  aF{  (5)  jp'<*  (23b) 


_  J  *2 

( s)p  cos  aFs  (s)p'ds 


(23c) 
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(24c) 


As  we  vary  <j>'  while  holding  s  constant  (e.g.  p',z'  constant),  R{<!>') is  an  even 


function:  R(-<f>')  =  R($')  ■  Therefore 


is  also  an  even  function  of  <f>' . 


Since  sin^'  is  an  odd  function  of  ^'the  integrand  in  Equation  (24a)  is  also  odd  giving 
Fs  ( s )  s  0.  We  are  left  with  only  Hy  in  (10b)  since  Hx  =  Hz  =  0 . 

Extending  this  to  any  position  r ,  we  have  for  Js  =  J,(s)t , 

H  =  (25a) 

where 

=  —  J  J,  (-s)|[yo'  sin  a  -  (z  -  z')  cos a]Fc  (5)  -psin  aFx  (s^p'ds  (25b) 

l7i 


(26a) 


(26b) 


E.  CENTERED  SPHERE  INTEGRATION  TEST  CASES 

A  centered  sphere  was  chosen  as  the  test  case  to  determine  the  accuracy  of  the 
integration  algorithm.  The  integration  algorithm  was  developed  using  the  trapezoidal  rule. 
Nine  cases  were  tested  with  sphere  sizes  ranging  from  a  radius  of  one  wavelength  to  ten 

wavelengths.  For  each  sphere  size  the  H  field  was  generated  at  three  different  locations: 
0.2,  1  and  3  wavelengths  from  the  outer  surface  of  the  sphere. 
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The  offset  distance  of  the  dipole  was  kept  constant  for  each  of  the  three  cases  associated 
with  the  three  different  sphere  sizes.  The  number  of  modes  used  was  determined  by  the 
following: 

N=integer(2  ka  +  2  ),  (27) 

where  k  is  the  wave  number  (2k  /  A)  and  a  is  the  sphere  radius.  In  each  of  the  nine 
cases  the  frequency  of  operation  was  300MHz.  The  RMS  error  was  determined  by 
comparing  the  H  field  determined  through  integration  to  the  H  field  exact  solution  found 
using  the  procedure  developed  in  Section  D.C.  Each  of  the  cases  was  tested  with 
increasingly  fine  segmentation  in  6  and  <f>  on  the  sphere  surface  until  the  RMS  error  was 
in  the  neighborhood  of  one  percent.  In  order  to  achieve  the  desired  integration  accuracy, 
the  number  of  integration  points  on  the  body  was  increased.  This  relates  to  the  spatial 
resolution  required  when  measuring  the  H  field  in  order  to  achieve  acceptable  results 
when  back  propagating  that  field.  As  the  number  of  integration  points  on  the  body 

increases,  the  measured  H  field  becomes  more  accurate.  The  more  accurate  the 
measured  field,  the  better  the  surface  current  on  the  sphere  can  be  resolved.  The  test 
cases  summarized  in  Table  1  are  chosen  based  solely  on  the  RMS  error  associated  with 
each  being  in  the  neighborhood  of  one  percent. 

The  results  shown  in  Table  1  indicate  the  RMS  error  for  source  integrated//^ 

observed  at  the  “field  distance”  radius.  The  RMS  error  is  a  function  of  the  number  of 
surface  points  at  which  the  surface  current  was  integrated.  Figures  9  through  12  show 
how  the  integration  results  converge  for  several  of  the  cases  shown  in  Table  1 .  For  the 
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CENTERED  SPHERE  INTEGRATION  RESULTS 

Case 

Sphere 

Radius 

<*) 

Field 

Distance 

(X) 

Segments 

e 

Segments 

Dipole 

Offset 

<*) 

Modes 

H* 

RMS  Error 

(%) 

la 

1 

1.2 

21 

21 

5 

14 

1.1490 

lb 

1 

2 

21 

21 

5 

14 

1.1620 

lc 

1 

4 

21 

21 

5 

14 

1.1510 

2a 

5 

5.2 

78 

78 

9 

64 

1.0310 

2b 

5 

6 

78 

78 

9 

64 

0.1220 

2c 

5 

8 

78 

78 

9 

64 

0.1250 

3a 

10 

10.2 

157 

157 

14 

128 

0.8268 

3b 

10 

11 

157 

157 

14 

128 

0.0456 

3c 

10 

13 

157 

157 

14 

128 

0.0459 

Table  1.  Centered  Sphere  Integration  Results 


cases  in  Table  1,  the  number  of  <f>  and  $  segments  is  based  on  the  radius  of  the  sphere. 
However,  the  number  of  <t>  segments  remains  constant  as  the  radius  of  each  ring  (see 
Figure  7)  decreases  near  both  poles  of  the  sphere.  This  finer  (f>  segmentation  gives 
increasingly  narrower  patches  of  surface  area  over  which  the  integration  is  performed. 
Consequently,  the  RMS  error  does  not  appear  to  decrease  significantly  as  the  number  of 
<f>  segments  is  increased.  This  is  a  result  of  the  increased  <f>  segmentation  that  is  “built  in” 
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as  the  radius  of  the  rings  decrease.  An  increase  in  RMS  error  is  observed  for  the 
decreased  field  distance  from  the  sphere  in  Cases  2a  and  3  a.  Since  the  field  distance  is 
only  0.21  from  the  sphere  surface  in  these  cases,  the  large  relative  contribution  and  rapid 
variation  of  inverse  distance  terms  in  the  integrand  of  Equation  (18)  for  surface  points 
closest  to  the  field  point  require  increasingly  fine  segmentation  for  accurate  numerical 
integration. 


Integration  Convergence:  Casela 
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Integration  Convergence:  Casel  b 


10 


Integration  Convergence:  Case2a 


70 
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Integration  Convergence:  Case2b 


90  90 

theta  segments 


phi  segments 


Figure  12.  Integration  Convergence:  Case  2b 


The  magnitude  and  phase  of  the  H  field  at  the  measured  distance  are  presented  in 

Figures  13  through  30.  These  figures  help  to  tell  the  complete  story  of  the  computed  H 
field.  In  each  figure  the  “exact”  quantity  computed  via  series  solution  is  depicted  by  a  line 
and  the  “measured”  quantity  is  depicted  by  dots.  Note  that  in  both  the  magnitude  and  the 
phase  plots  the  measured  quantity  closely  matches  the  exact  quantity.  These  results 

indicate  that  the  integration  produces  an  H  field  with  minimal  error.  This  knowledge  will 
help  isolate  the  source  of  error  that  may  present  itself  due  to  the  back-propagation 
process,  which  will  be  described  in  Chapter  III. 
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Phase  Angle  (radians) 


Theta  (degrees) 


Figure  13.  Magnitude  Comparison:  Case  la 


Hexact  (line)  &  Hinteg  (dots):  Casela 
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Hexact  (line)  &  Hinteg  (dots):  Caselc;  RMS  Error=  1.151% 


Hexact  (line)  &  Hinteg  (dots):  Caselc 
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Phase  Angle  (radians) 


Hexact  (line)  &  Hinteg  (dots):  Case2a;  RMS  Error=  1.031% 


Hexact  (line)  &  Hinteg  (dots):  Case2a 
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Phase  Angle  (radians) 


Hexact  (line)  &  Hinteg  (dots):  Case2b;  RMS  Error=  0.122% 


Hexact  (line)  &  Hinteg  (dots):  Case2b 
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Phase  Angle  (radians) 


Hexact  (line)  &  Hinteg  (dots):  Case3b;  RMS  Error=  0.04562% 


Hexact  (line)  &  Hinteg  (dots):  Case3b 
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Figure  30.  Phase  Comparison:  Case  3c 


m.  BACK  PROPAGATION  OF  SCATTERED  FIELD  (MEASURED  FIELD) 


A.  BACKGROUND 

A  finite  element  method  is  used  to  back-propagate  the  measured  field  to  the 
surface  of  the  sphere.  The  measured  fields  that  are  back-propagated  are  those  determined 
in  Chapter  II  and  summarized  in  Table  1.  Once  the  measured  field  is  back-propagated,  the 
surface  current  is  determined  using  Equation  (9b).  The  surface  current  due  to  back- 
propagation  (Je^  )  is  compared  to  the  original  surface  current  ( Je^ )  due  to  the 
elemental  dipole  which  is  computed  using  a  spherical  harmonic  series,  per  Section  II.C. 
The  following  section  describes  the  finite  element  formulation. 

B.  FINITE  ELEMENT  FORMULATION  FOR  AXISYMMETRIC  FIELDS 

For  the  special  case  being  considered,  with  no  <t>  variation  in  materials  or  sources, 
and  Jf  =  0,  the  fields  are  derived  in  Section  H.B.,  results  of  which  are  repeated  here  for 


convienence, 

H(p,z)  =  Ht(p,z)j  (28a) 

E(p,z)  =  Et(p,z)p+E,(p,  z)z  (28b) 

Maxwell’s  curl  equations  are:  (with  <y  =  2 jrf) 

V  x  7/  =  jasE  +  J  (29a) 

Vx£  =  -  ja/iH  (29b) 


37 


Using  Equation  (29a)  for  assumed  fields  in  Equation  (28)  gives 


•  zr  ~im*  1  T 


(30a) 


1 1  a  i 

]coE’=~e~p'7p{pH>)~~eJ- 


(30b) 


These  are  the  generating  equations  for  E  in  terms  of  H+  and  known  source 


current  J .  In  this  case  Equation  (29b)  simplifies  to: 


-  ~  cE  cE2  , 

j<o(VxE)-<j>  =jco  =0) 


Substituting  the  E  fields  from  Equation  (30)  into  Equation  (3 1)  gives  the  partial 
differential  equation  (PDE)  satisfied  by  Hy. 


<?[  d  ,  d\\&lA  2  rj  &  A  j\  d  A  n 

dp  ep  dp +  dz &  J ^ + ' ®  fJ1*  ~  dp jz)  dz^s  Jp) 


where  s  =  s(p,  z)  is  allowed  to  be  to  be  inhomogeneous. 

The  PDE  in  Equation  (32)  will  be  numerically  solved  using  the  finite  element 
method  (FEM),  where  is  approximated  with  electrically  small  triangular  element 

regions  in  the  (p,z)  plane.  [12] 

A  variational  Euler-Lagrange  algorithm  is  used  to  implement  the  finite  element 
method.  This  algorithm  replaces  the  solution  of  the  PDE  in  Equation  (32)  with  finding  the 
stationary  point  in  complex  function  space  of  a  quadratic  (in  H , )  functional.  We  first 

expand  the  field  using  pyramidal  basis  functions,  un(p,z),  defined  for  each  node  in  the 


mesh  and  having  support  of  the  surrounding  triangular  elements, 


38 


(33) 


N 

H+(p,z)  =  2>„w„(p,r) 

n~\ 

and  then  substitute  Equation  (33)  into  the  functional.  The  stationary  solution  is  found  by 
differentiating  with  respect  to  the  unknown  values  of  ,  located  at  the  m-th  nodes,  to 

yield  the  linear  system 

2>„  jj  ^  V(pwm)  •  V(pw„)  -  k2pumu„dpdz  =  0  (34) 

n=l  overlapping  elements  • 

Defining  the  double  integrals  indexed  on  m,n  as  T (m,ri) ,  the  linear  system  in  Equation 
(34)  can  be  rewritten  as, 

Z  hj(m,n)  =  ~  ^h„T(m,n)  (35) 

n  inside  boundary  n  on  boundary 


->na(w) 

Hb(n) 

0" 

"  K  ~ 

i 

na(m) 

A 

° 

K 

J1 Nm . 

= 

B 

Nm  xN*  Square  Am# 

( very  sparse  ) 

Nm*Ns 

{sparse) 

boundary 

values 


(36) 


The  m-th  row  of  A  and  B  is  formed  by  Equation  (35)  with  the  n-th  term 
corresponding  to  a  column  of  A  if  hn  is  unknown.  The  appropriate  column  of  B  is  filled 
if  hn  is  a  boundary  node.  There  is  no  contribution  to  A  or  B  for  nodes  on  the  z-axis 
(p  =  0)  since  H+(p  =  0)  =  0  and  those  nodal  values  of  hn  are  known.  [13] 


C.  FINITE  ELEMENT  METHOD  LIMITATIONS 

The  finite  element  method  was  successful  for  back-propagation  of  fields  measured 
at  distances  less  than  0.25  wavelengths  from  the  surface  of  the  sphere.  When  the  field  was 
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measured  at  distances  greater  than  0.25  wavelengths,  errors  due  to  back-propagation 
increased  significantly.  The  same  cases  for  which  the  measured  field  was  generated  were 
used  for  back-propagation.  These  results  are  summarized  in  Table  2.  The  program 
FEM2.M  was  used  to  determine  the  surface  current  due  to  back-propagation.  It  can  be 
found  in  the  Appendix. 


BACK-PROPAGATION  RESULTS:  FINITE  ELEMENT  METHOD 

Case 

Sphere 

Field 

Modes 

H+ 

Radius 

Distance 

RMS  Error 

RMS  Error 

(to 

'(to 

(%> 

(%) 

FEMla 

1 

1.2 

14 

0.0590 

0.9976 

FEMlb 

1 

2 

24 

0.0338 

30.09 

FEMlc 

1 

4 

50 

0.0147 

25.87 

FEM2a 

5 

5.2 

64 

0.0043 

0.7406 

FEM2b 

5 

6 

74 

0.0052 

57.18 

FEM2d 

5 

5.26 

64 

0.0042 

662.9 

FEM3a 

10 

10.2 

128 

0.0017 

0.7333 

FEM3d 

10 

10.26 

128 

0.0063 

5.547 

Table  2.  Centered  Sphere  Back-Propagation  Results:  Finite  Element  Method 
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The  magnitude  and  phase  of  the  surface  current  on  the  PEC  sphere,  as  determined 
through  the  finite  element  method,  is  shown  in  Figures  3 1  through  42.  These  figures  help 
to  tell  the  complete  story  of  the  accuracy  in  back  propagating  to  obtain  the  current.  In 
each  figure  the  “exact”  quantity  is  depicted  by  a  solid  line  and  the  calculated  quantity  is 
depicted  by  dots.  Note  that  in  both  the  magnitude  and  phase  plots  the  calculated  quantity 

closely  matches  the  exact  quantity  for  the  cases  in  which  the  H  field  is  measured  less  than 
0.25  wavelengths  away  from  the  surface  of  the  sphere.  In  the  cases  where  the  distance  is 
one  wavelength,  the  error  is  immense.  Cases  FEM2d  and  FEM3d  have  been  added  to  the 
original  nine  test  cases  in  order  to  demonstrate  what  happens  when  the  measurement 
distance  is  increased  to  0.26  wavelengths  from  the  surface  of  the  sphere.  The  following 
section  will  investigate  the  cause  of  this  error. 
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Phase  Angle  (radians) 


Jexact  (line)  &  Jback  (dots):  FEMIa;  RMS  Error=  0.9976% 


Jexact  (line)  &  Jback  (dots):  FEMIa 
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Figure  34.  Phase  Comparison:  FEMlb 
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Jexact  (line)  &  Jback  (dots):  FEM2a;  RMS  Error-  0.7406% 


Jexact  (line)  &  Jback  (dots):  FEM2a 


Theta  (degrees) 

Figure  38.  Phase  Comparison:  FEM2a 


Phase  Angle  (radians) 


Jexact  (line)  &  Jback  (dots):  FEM2b;  RMS  Error=  57.18% 


Jexact  (line)  &  Jback  (dots):  FEM2b 
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Jexact  (line)  &  Jback  (dots):  FEM2d;  RMS  Error=  662.9% 


Jexact  (line)  &  Jback  (dots):  FEM2d 
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Phase  Angle  (radians) 


Jexact  (line)  &  Jback  (dots):  FEM3a;  RMS  Error=  0.7333% 


Theta  (degrees) 


Figure  43.  Magnitude  Comparison:  FEM3 a 


Jexact  (line)  &  Jback  (dots):  FEM3a 
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Figure  46.  Phase  Comparison:  FEM3d 
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When  back-propagating  measurements  are  made  at  distances  greater  than  0.25 
wavelengths  from  the  surface  of  the  sphere  using  the  FEM  approach,  the  RMS  error  of  the 
predicted  surface  current  tends  to  increase  dramatically  relative  to  the  situation  of  smaller 
distances.  This  appears  to  result  from  resonant  modes  injecting  themselves  into  the 
solution  in  varying  degrees.  Figure  47  shows  the  single  mode  RMS  error  for  a  centered 
PEC  sphere  used  in  case  FEM  la.  Note  that  the  RMS  error  associated  with  each  mode  is 
relatively  low;  all  but  two  are  below  one  percent.  The  results  shown  Table  2  indicate  that 
back-propagation  gave  a  solution  with  less  than  one  percent  RMS  error.  Figure  48  shows 
the  single  mode  RMS  error  for  case  FEMlb.  In  this  case  the  measured  field  was  one 
wavelength  from  the  surface  of  the  sphere.  Note  that  the  lower  order  modes  generally 
have  RMS  error  in  the  neighborhood  of  five  percent  or  less,  except  for  mode  n  =  5  which 
is  nearly  fifty  percent  RMS  error.  This  error  “spike”  indicates  a  resonant  mode  that 
significantly  alters  the  solution.  As  can  be  seen  in  Figures  48  through  53,  the  resonant 
modes  seem  to  appear  in  only  the  cases  where  the  measurement  distance  is  greater  than 
0.25  wavelengths  from  the  surface  of  the  sphere. 
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Percent 


Single  Mode  RMS  Error  for  Centered  PEC  Sphere:  FEMIa 


Figure  47.  Single  Mode  RMS  Error:  Case  FEMIa 


Single  Mode  RMS  Error  for  Centered  PEC  Sphere:  FEMIb 


Figure  48.  Single  Mode  RMS  Error:  Case  FEMIb 


51 


Mode  Index  "n" 

Figure  49.  Single  Mode  RMS  Error:  Case  FEM2a 


Single  Mode  RMS  Error  for  Centered  PEC  Sphere:  FEM2b 


Single  Mode  RMS  Error  for  Centered  PEC  Sphere:  FEM3a 


Single  Mode  RMS  Error  for  Centered  PEC  Sphere:  FEM3d 


Figure  53.  Single  Mode  RMS  Error:  CaseFEM3d 

Another  way  to  view  this  phenomenon  is  by  looking  at  the  product  solution  parts 
of  H^(r,  6)  =  g(r )  •  h(0)  for  a  given  mode,  n.  The  radial  variation,  g(r) ,  to  be  shown  in 

the  upper  subplots,  indicates  the  contribution  at  the  PEC  sphere  surface  due  to  a  constant 
contribution  of  g(b)  =  1  on  the  outer  boundary.  Figures  54  and  55  show  the  radial 
variation  for  modes  n  =  4  and  n  =  5,  respectively,  for  the  case  FEMlb.  Note  that  mode  4 
provides  relatively  small  contribution  at  the  PEC  sphere  surface.  Mode  5,  however, 
indicates  an  extremely  large  contribution  at  the  PEC  sphere  surface.  Two  additional  cases 
have  been  added  here  to  further  illustrate  this  observation.  Cases  FEM2d  and  FEM3d 
shown  in  Figures  56  and  57,  respectively,  are  both  cases  in  which  the  measured  field  is 
0.26  wavelengths  from  the  surface  of  the  sphere.  Note  the  large  contribution  at  the 
surface  of  the  sphere  due  to  the  resonant  modes  captured  in  each  of  these  figures. 
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Radial  Variation 


Product  Solution  Parts  of  Hphi(r, theta):  FEM2d;  n=8 


Product  Solution  Parts  of  Hphi(r, theta):  FEM3d;  n=16 


Figure  57.  Product  Solution  Parts:  Case  FEM3b,  n  =  16 
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Additional  insight  is  provided  by  viewing  the  resonant  modes  that  appear  when  the 
field -is  measured  at  varying  distances  from  a  fixed  size  sphere.  Figure  58  considers  modes 
n  =  1  to  100.  Each  line  in  Figure  58  displays  the  contribution  of  g(a)  from  a  single  mode 
at  the  sphere  surface  as  a  result  of  a  constant  contribution  of  g(b )  =  1  at  the  boundary 
distance.  The  “best  radius,”  b  =  5.641  ^  in  this  case,  provides  the  minimum  g(a)  from 
all  modes  applied  as  g(b )  =  1  .  From  the  perspective  of  minimal  resonance  effects,  this  is 
the  best  location  to  measure  the  field. 


b  (wavelengths) 

Figure  58.  Resonant  Effects  at  Various  Outer  Boundary  Locations 

D.  TRANSFER  MATRIX  FOR  CENTERED  SPHERE 

As  can  be  seen  in  Table  2  the  results  due  to  back-propagation  are  poor  for 
measurements  made  greater  than  approximately  0.25  wavelengths  from  the  surface  of  the 
sphere.  A  transfer  matrix  using  spherical  harmonics  for  a  centered  sphere  was  developed 
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to  accurately  back  propagate  the  field  on  the  outer  boundary  to  that  on  the  sphere, 
without  the  bothersome  injection  of  resonant  modes  as  appear  with  the  FEM  solution. 

This  method  will  give  some  indication  as  to  the  spatial  resolution  required  when  measuring 
the  field  so  that  the  source  current  on  the  sphere  can  be  accurately  determined. 

Figure  58  shows  a  meridian  plane  for  a  centered  sphere  of  radius  a  enclosed  by  an 

outer  boundary  arc  of  radius  b .  For  any  radius  r ,  where  a  <r  <b ,  the  total  H  field  can 
be  determined  by:  [14] 


N 


Ht(r,S)  £  2>, 


«=  1 


w+b 


kr 


kr 


P„'(cosff) 


(37) 


where: 


forces  Eg  =0  for  each  mode  on  the  surface  of  the  sphere. 


(38) 


A  specified  f(9)  =  H+(by0)  will  provide  a  set  {a„} .  The  an  coefficients  are 

determined  by  moment  matching,  using  orthogonality  of  Legendre  functions.  This 
technique  is  detailed  in  [14  ].  The  transfer  matrix  Hs is  developed  such  that 


Hf(a,Oq)  =  Hs(q,p)- Hf(b,0p) 


(39) 
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Centered  Sphere  radius  =  2;  Arc  radius  =  5 


Figure  59.  Meridian  Plane  of  Centered  Sphere 

where  Ht{b,6p )  is  a  column  array  of  the  0p  sampled  field  at  radius  b ,  as  determined  in 
Chapter  II.  The  array  product  results  in  the  column  array  of  at  radius  a  which 

corresponds  to  the  surface  of  the  sphere.  From  this  result  the  surface  current  is  then 
determined  using  Equation  (9b).  The  test  cases  using  the  transfer  matrix  were  generated 
using  the  program  HSPHERE1.M  which  can  be  found  in  the  Appendix.  Table  3  tabulates 
the  results  and  compares  them  to  those  obtained  using  the  finite  element  method.  It 
should  be  noted  that  the  transfer  matrix  works  only  for  a  centered  sphere  and  can  not  be 
applied  to  any  arbitrary  axisymmetric  body. 
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SURFACE  CURRENT  RESULTS  USING  TRANSFER  MATRIX 


Case 

Sphere 

Radius 

0) 

Field 

Distance 

<W 

Modes 

e 

Segments 

H* 

RMS  Error 

(%) 

J0 

RMS  Error 

(%) 

la 

1 

1.2 

14 

94 

0.0590 

0.8085 

lb 

1 

2 

24 

125 

0.0338 

0.6271 

lc 

1 

4 

50 

188 

0.0147 

0.7349 

2a 

5 

5.2 

64 

408 

0.0043 

0.7946 

2b 

5 

6 

74 

376 

0.0052 

0.9453 

2d 

5 

5.26 

64 

408 

0.0042 

0.7764 

3a 

10 

10.2 

128 

801 

0.0017 

0.7885 

3d 

10 

10.26 

128 

801 

0.0063 

0.7827 

Table  3 .  Surface  Current  Results  Using  Transfer  Matrix 


The  magnitude  and  phase  of  the  surface  current  on  the  PEC  sphere,  as  determined 
by  the  transfer  array,  is  shown  in  Figures  60  through  77.  Once  again  these  figures  help  to 
tell  the  complete  story  of  the  surface  current.  In  each  figure  the  exact  quantity  is  depicted 
by  a  solid  line  and  the  calculated  quantity  is  depicted  by  dots.  Note  that  in  both  the 
magnitude  and  phase  plots  the  calculated  quantity  closely  matches  the  exact  quantity. 
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Phase  Angle  (radians) 


Jexact  (line)  &  Jback  (dots):  Casela;  RMS  Error=  0.8085% 


Jexact  (line)  &  Jback  (dots):  Casela 
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Jexact  (line)  &  Jback  (dots):  Caselc;  RMS  Error=  0.7349% 


Jexact  (line)  &  Jback  (dots):  Caselc 


Phase  Angle  (radians) 


Jexact  (line)  &  Jback  (dots): 

0.1 1 - t - 
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Jexact  (line)  &  Jback  (dots):  Case2a 
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Phase  Angle  (radians) 
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Figure  71.  Phase  Comparison:  Case2c 
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Jexact  (line)  &  Jback  (dots):  Case3a;  RMS  Error=  0.7885% 
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Figure  72.  Magnitude  Comparison:  Case3a 


Jexact  (line)  &  Jback  (dots):  Case3a 
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Figure  75.  Phase  Comparison:  Case3b 


200 


68 


Jexact  (line)  &  Jback  (dots):  Case3c;  RMS  Error=  0.7846% 


Theta  (degrees) 


Figure  76.  Magnitude  Comparison:  Case3c 


Jexact  (line)  &  Jback  (dots):  Case3c 


Theta  (degrees) 

Figure  77.  Phase  Comparison:  Case3c 
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IV.  ANALYSIS  OF  RESULTS 


A.  THE  INTEGRATION  ALGORITHM 

Since  the  “measured”  fields  to  be  back-propagated  were  to  be  determined  through 
integration,  the  first  part  of  this  thesis  was  dedicated  to  generating  accurate  fields  through 
numerical  integration.  As  is  shown  in  Table  1,  properly  performed  numerical  integration 
provides  highly  accurate  measured  fields.  Figures  9  through  12  provide  several  examples 
of  how  quickly  the  integration  converges.  In  each  case  the  percent  RMS  error  is 
determined  essentially  by  the  number  of  points  on  the  sphere  surface  over  which  the 
integration  takes  place.  For  the  cases  included  in  this  thesis  the  number  of  points  was 
chosen  such  that  the  error  due  to  back-propagation,  as  determined  by  the  transfer  matrix, 
would  be  kept  below  one  percent. 

Figures  13  through  30  show  the  magnitude  and  phase  of  the  measured  field 
determined  through  numerical  integration.  For  each  figure  the  measured  field  is  over-laid 
on  the  exact  field.  Visual  inspection  indicates  the  magnitudes  are  nearly  the  same  over  the 
entire  range  of  theta.  RMS  error  computations  indicate  a  difference  of  only  a  few 
hundredths  of  one  percent.  Visual  inspection  of  the  phase  comparison  gives  similar  results 
-  both  the  computed  and  the  exact  are  nearly  the  same. 

As  the  number  of  surface  points  increase  for  a  constant  sphere  size,  the  accuracy 

of  the  measured  H  field  increases.  This  result  takes  the  form  of  decreasing  RMS  error  as 
shown  in  the  integration  convergence  of  Figures  9  through  12.  As  the  measurement 
distance  increases  for  a  constant  sphere  size,  the  number  of  points  required  for  accurate 
integration  results  decreases.  This  can  be  seen  by  comparing  Cases  2a  and  2b  shown  in 
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Figures  1 1  and  12.  For  the  larger  measurement  distance  of  Case  2b  the  RMS  error  is 
significantly  smaller  using  the  same  number  of  points  as  Case2a. 


B.  FINITE  ELEMENT  METHOD  PROBLEMS 

The  finite  element  method  was  used  successfully  to  determine  the  source  current 
on  the  PEC  sphere  when  the  field  was  measured  at  distances  less  than  0.25  wavelengths 
from  the  sphere  surface.  Table  2  summarizes  the  test  cases  used  in  this  thesis. 

The  large  errors  that  arise  when  the  field  to  be  back-propagated  was  measured  at 
distances  greater  than  0.25  wavelengths  from  the  sphere  surface  appear  to  be  caused  by 
resonant  modes.  These  modes  are  shown  for  several  cases  in  Figures  48,  50,  52  and  53. 
When  back-propagated,  these  modes  provide  overwhelmingly  disproportionate 
contributions  at  the  sphere  surface.  These  contributions  are  shown  for  a  few 
representative  cases  in  Figures  55,  56,  and  57.  It  appears  these  resonant  modes  come 
about  as  a  result  of  the  boundary  conditions  placed  on  the  outer  boundary,  the  location  at 
which  the  field  is  measured.  As  can  be  seen  in  Figures  48,  50,  51  and  53,  these  resonant 
modes  appear  as  sharp  spikes  in  the  plot. 

C.  TRANSFER  MATRIX 

The  transfer  matrix  that  was  used  to  determine  the  source  current  on  the  body 
gave  results  with  less  than  one  percent  RMS  error.  The  low  error  is  a  result  of  the  number 
of  surface  points  used  to  determine  the  measured  field,  and  the  number  of  points  at  which 
the  field  was  measured  at  the  outer  boundary.  The  relationship  between  number  of 
surface  points  and  numerical  integration  accuracy  was  discussed  in  Section  A.  As  can  be 
seen  in  Equation  (18),  the  measured  field  varies  inversely  with  distance  from  segment  to 
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field  point.  This  inverse  distance  variation  is  the  main  point  of  concern  when  R  is  small  at 
the  closest  point  of  approach.  There  is  a  large  relative  contribution  to  the  integration 
from  this  surface  region  and  the  variation  in  the  field  is  rapid.  The  number  of  segments  on 
the  sphere  surface  in  this  near-region  must  be  increased  to  ensure  the  accuracy  of 
numerical  integration. 
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V.  CONCLUSIONS 


The  objective  of  this  research  was  to  investigate  the  viability  of  back-propagating 
electromagnetic  field  measurements  to  an  axisymmetric  body  in  order  to  determine  the 
source  distribution  for  radiated  power  from  the  body.  The  method  for  achieving  this 
objective  was  straightforward.  First,  the  “measured  fields”  were  simulated  through 
numerical  integration.  Second,  these  fields  were  input  into  a  finite  element  algorithm  to 
determine  the  source  current  on  a  PEC  axisymmetric  body.  The  original  thrust  of  the 
thesis  was  to  investigate  the  ranges  and  resolution  at  which  the  fields  could  be  measured  in 
order  to  provide  a  usable  level  of  accuracy  in  determining  the  source  distribution  on  the 
body.  The  test  cases  were  to  include  various  “arbitrary”  axisymmetric  shapes  to  include 
cones,  offset  cones,  offset  spheres,  and  various  cone-cylinder-sphere  combinations.  The 
centered  sphere  was  to  be  used  for  very  basic  test  purposes  only. 

Along  the  way  to  the  objective  it  was  found  that  the  finite  element  method 
provided  inaccurate  results  in  the  centered  sphere  tests  for  all  but  very  limited  cases.  This 
discovery  prompted  a  redirection  of  efforts  to  investigate  the  sources  of  error. 
Consequently,  progress  towards  the  original  goals  was  sidetracked;  but,  that  is  often  the 
nature  of  true  research  into  the  unknown. 

The  integration  program  developed  to  determine  the  measured  fields  was  found  to 
converge  to  an  accurate  solution.  Although  the  shape  tested  was  a  sphere,  the  theory 
developed  in  Chapter  II  using  stacked  circular  rings  can  be  easily  extended  to  arbitrary 
axisymmetric  shapes. 

The  finite  element  method  was  selected  as  a  means  to  determine  the  source  current 
on  the  sphere  through  back-propagation.  The  errors  associated  with  this  method  for  the 
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centered  sphere  were  unexpected.  Consequently,  solving  these  associated  resonance 
problems  was  not  intended  to  be  the  focus  of  this  thesis.  Further  testing  of  the  FEM 
solution  will  be  done  as  part  of  extensions  to  this  effort  in  follow-on  thesis  efforts  and 
techniques  for  mitigating  resonance  effects  will  be  explored. 
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APPENDIX  -  MATLAB  CODES 


%thesis2.M  by  Dan  Wawrzyniak  5/10/97 
%  updated  5/14/97 
%  Combines  all  programs  into  one 
%  Includes  checking  integration  error  first 
%  Allows  user  to  save  data  in  a  .mat  file  for  future  use 
clear  all; 

case  =  input('What  case  are  you  running?  ’/s'); 
fO  =  input('Enter  frequency  in  MHz: '); 
lambda  =  3e2/f0; 

disp(['Wavelength  =  ',num2str(lambda),'  meters']); 

al  =  input('Enter  metal  sphere  radius  in  terms  of  wavelength:  ');a  =  al  *lambda; 
rl  =  input(Enter  outer  boundary  in  terms  of  wavelength:  ');r  =  rl  “"lambda; 
eta  =  120*pi;  %  free-space  Z_0 

Idl  =  1;  %  dipole  moment 

B  =  2*pi/lambda;  %  wavenumber 

zOl  =  input(Enter  z-position  of  elemental  dipole  in  terms  of  wavelength: '); 
zO  =  z01*lambda; 

Nthl  =  input(Enter  number  of  theta  points: ');  Nth  =  Nthl-1; 
ppwl  =  input(Enter  number  points  per  wavelength: '); 
dth=pi/Nth;  th=(0:dth:pi)';  thdeg=180*th/pi; 
sth=sin(th);  cth=cos(th); 

%  computing  rd,  sin(thd),  and  cos(thd)  for  sphere  surface  (r=a) 
rd=sqrt(a*a+zO*zO-2*a*zO*cth);  sthd=a*sth./rd;  cthd=(a*cth-zO)./rd; 

%  computing  field  components  in  dipole  centered  coordinates 
Fth=Idl*exp(-j*B*rd)./(4*pi*rd); 

HpHFth.  *sthd.  *(j  *B  +  1  ,/rd); 

Etd=eta*Fth.  *sthd.  *(j*B  +  l./rd  -  j./(B*rd.*rd)); 

Erd=2*eta*Fth.*cthd.*(l./rd  -  j./(B*rd.*rd»; 

%  translating  incident  field  spherical  components 
cdth=cth.*cthd  +  sth.*sthd;  sdth=sth.*cthd  -  cth.*sthd; 

Eti=Etd.*cdth  -  Erd.*sdth;  %  note  that  Hpi=Hpd 
%  Computing  spherical  harmonic  coefficients  for  scattered  field 
%  Ets  expansion  which  minimizes  the  total  sphere  surface 
%  tangential  field  =  Eti+Ets 
Nmax  =  fix(2*B*r); 
disp([Nmax  =  ',num2str(Nmax)]); 

N  =  input('Enter  number  of  spherical  harmonics: '); 

Sn=-dth*Eti.*sth;  Pnl=legpol2(cth,N);  In=(Pnl.')'"Sn; 
n=(l  :N)';  cn=(2*n+l)./(2*n.*(n+l»; 

Ba=B*a;  [Hn,DHn]=shan3(Ba,N); 

Hn  =  Hn.';  DHn  =  DHn.'; 

Dn=DHn/Ba;  an=(cn.*In)./(j*eta*Dn); 

Ets=j*eta*Pnl  *(an.  *DHn)/Ba;  Ett=Eti+Ets;  %  ideally  Ett  -->  0 
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Hps=Pnl*(an.*Hn)/Ba;  Hpt=Hpi+Hps;  %  Jtheta=H_pt 

Jtheta  =  -Hpt; 

Br=B*r;  [Hn,DHn]=shan3(Br,N); 

Hn  =  Hn.';  DHn  =  DHn.'; 

Ets=j  *eta*Pnl  *(an.  *DHn)/Br;  Hps=Pn  1  *(an.  *Hn)/Br; 
numring  =  Nth;  %  number  of  rings 
inc  =  inputCEnter  number  of  face  segments: '); 
ppoints  =  input('Enter  phipoints  -  wavelength  factor: '); 

%  Field  Point  Inputs 

%  Input  the  Field  Point  in  terms  of  theta, 
radius  =  r;  %  radius  of  field  point 
numFP  =  length(Hps); 
thetaFP  =  linspace(0, pi, numFP); 
thetaFPd  =  180*thetaFP/pi; 
rho  =  radius.  *sin(thetaFP); 
z  =  radius.  *cos(thetaFP); 

%  Jt  is  the  tangential  surface  current  at  the  nodes. 

%  The  nodes  are  defined  by  rhop  &  zp. 

Jt  =  Jtheta; 

%  Determine  the  thetas  for  the  sphere 

dtheta  =  pi/numring;  %  gives  the  delta  theta 

theta  =  0:dtheta:pi;  %  gives  numring  +1  values  of  theta. 

thetad  =  theta.  *  1 80/pi;  %  theta  in  degrees 

%  Determine  (rhoprime,zprime)  pairs  on  the  sphere 
rhop  =  a*sin(theta); 
zp  =  a*cos(theta); 

%  Find  the  length  "s"  of  each  face. 

s  =  sqrt(rhop(2)A2  +  (zp(l)-zp(2))A2); 

%  Find  the  angle  associated  with  each  face 

alpha  =  acos(rhop(2)/s)  +  theta(l  :numring);  %  radians 
alphad  =  alpha*  1 80/pi;  %  degrees 

talpha  =  tan(alpha); 
alpha  1  =  ones(inc+ 1 , 1  )*alpha; 
alpha  =  reshape(alphal,numring*(inc+l),l); 
dzf  =  (zp(l:numring)-zp(2:numring+l))/inc; 
phipoints  =  round(pi.*a*ppoints/lambda)+l; 
dphi  =  pi./phipoints; 
phi  =  0:dphi:pi; 
cp  =  cos(phi); 

for  n  =  1  :numring;  %  loops  through  rings 
zf(:,n)  =  (zp(n) :  -dzf(n) :  zp(n+ 1 ))'; 
rhof(:,n)  =  (linspace(rhop(n),rhop(n+l),inc+l))'; 

JT(:,n)  =  Jt(n)*(zf(l:inc+l,n)-zp(n+l))/(zp(n)-zp(n+l))+... 

Jt(n+ 1  )*(zf(  1  :inc+ 1  ,n)-zp(n))/(zp(n+  l)-zp(n)); 
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%  Reshape  the  matrices  to  column  vectors 
zf  =  reshape(zf,numring*(inc+ 1 ),  1 ); 
rhof  =  reshape(rhof,numring*(inc+l),  1); 

JT  =  reshape(JT,numring*(inc+ 1 ),  1 ); 
for  FP  =  1  :numFP;  %  loops  through  field  points 
R1  =  2  *  rho(FP)  *  rhof*  cp ; 

R2  =  (rho(FP)A2+rhof.A2+(z(FP)-zf).A2)*ones(l,length(cp)); 

R  =  sqrt(R2-Rl); 

[Fes, FIs]  =  fcsfls2(cp,dphi,R); 

Ha  =  rhof.*sin(alpha); 

Hb  =  (z(FP)  -  zf).*cos(alpha); 

He  =  (Ha-Hb).*Fcs; 

Hd  =  rho(FP).*sin(alpha).*Fls; 

H  =  (Hc-Hd).*JT; 

H  =  reshape(H,inc+l,numring); 

HI  =  a*((l/(4*pi))*(trapz(H))); 

HI  =  dzf.*Hl; 

HH(l:numring,FP)  =  HI.'; 
end 

Hphi  =  sum(HH); 

%%  RMS  Error  Calculations 

E  =  100*sqrt((Hps  -  Hphi.')'*(Hps  -  Hphi.'))./sqrt(Hps'*Hps); 
figure 

plot(thdeg,abs(Hphi),'ko',thdeg,abs(Hps),'k') 

title(['Hps  (exact:  integrated:  o):  ',eval('case'),';  RMS  Error=  ',num2str(E),'%']) 

xlabel('Theta  (degrees)'); 

ylabel('|Hscattered|'); 

figure 

plot(thdeg,real(Hphi),'k*',thdeg,imag(Hphi),'k+,,thdeg,real(Hps),... 

'k',thdeg,imag(Hps),'ko'); 

title([Hexact(real:  imag:  o)  &  Hint(real:  *,  imag:  +):  ',eval('case')]) 
xlabel('Theta  (degrees)') 
ylabel('Hps  (Amps/meter)') 

%text(. 55,. 95,  ['Frequency=  ',num2str(fD),  'Mhz'], 'color', . . . 

%  'g','units','normalized'); 

%text(. 5 5,. 90,['Wavelength=  ’,num2str(lambda), 'meters’], 'color', . . . 

%  'g', 'units', 'normalized'); 

%text(.55,. 85, [Dipole  Offset=  ',num2str(z01),'*lambda'],'color',... 

%  'g', 'units', 'normalized'); 

%text(.55,. 80, [Radius  of  Sphere=  ',num2str(al),'*lambda'], 'color',... 

%  'g', 'units', 'normalized'); 

%text(.55,.75, [Radius  of  Field  Point=  ',num2str(rl),'*lambda'],'color',... 

%  'g','units','normalized'); 

%text(.55,.70,['Number  Rings=  ',num2str(numring)], 'color',... 

%  'g', 'units', 'normalized'); 
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%text(.55,.65,['Face  Segments= num2str(inc)], 'color',. . . 

%  'gVunitsVnormalized'); 

%text(. 55,  .60,  ['Phi-Wavelength  Factor=  ',num2str(ppoints)], 'color', . . . 

%  'g','units','normalized'); 

%text(.55,. 55, ['Field  Points=  ',num2str(numFP)],'color',... 

%  'g','units','normalized'); 

%text(.55,.50, [Field  Points  per  wl=  ',num2str(ppwl)], 'color',... 

%  'g','units','normalized'); 

%text(.55, .45, ['Number  modes=  ',num2str(N)],'color',... 

%  'g','units','normalized'); 

%text(.55,.40,[rRMS  Error=  ',num2str(E),’  %'], 'color',... 

%  'r','units','normalized'); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Computes  Hphi  incident  at  the  field  point  radius  due  to  the  dipole. 

%  This  result  is  added  to  the  Hphi  scattered  computed  using  Dan's 
%  integration  program 
dth=pi/Nth;  th=(0:dth:pi)';  thdeg=180*th/pi; 
sth=sin(th);  cth=cos(th); 

%  computing  rd,  sin(thd),  and  cos(thd)  for  sphere  surface  (r=a) 
rd=sqrt(r*r+z0*z0-2*r*z0*cth);  sthd=r*sth./rd;  cthd=(r*cth-zO)./rd; 

%  computing  field  components  in  dipole  centered  coordinates 
Fth=Idl*exp(-j*B*rd)./(4*pi*rd); 

Hpi=Fth.*sthd.*(j*B  +  l./rd); 

Hphi  =  Hphi.';  %  make  Hphi  a  column  vector 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  This  program  loads  the  transfer  function  Hs  from  femdat.mat 
%  and  created  by  FEM2.M  and  dots  it  with  Hphi  stored  in 
%  source.mat  from  Dan's  integration  program  SOURCE.M. 

% 

%  The  result  is  the  current  on  the  surface  of  the  sphere 
%  arrived  at  through  the  back  propagation  of  the  scattered 
%  H  field. 

% 

%  The  result  (Jback)  is  compared  to  the  surface  current  used  in 
%  Dan's  program  (Jtheta,  stored  in  jtheta.mat). 

% 

%  The  surface  current  is  generated  by  DIPSPHR2.M 
% 

load  femdat.mat 

%  Add  Hphi  scattered  from  sphere  and  Hpi  incident  from  dipole 
Hphitot  =  -(Hphi  +  Hpi); 

%  Determine  Surface  current  due  to  back-propagation  of  Hphi 
%  Hhpi  is  UNFILTERED 
Jback  =  Hs*  Hphitot; 

%  Determine  Surface  current  due  to  back-propagation  of  Hphi 
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%  RMS  Error  Calculations 

E  =  100*sqrt((Jtheta  -  Jback)'*(Jtheta  -  Jback))./sqrt(Jtheta'*Jtheta); 
figure 

plot(thdeg,abs(Jtheta),'k',thdeg,abs(Jback),'ko') 
title(['Jexact  (-)  &  Jback  (o):  ',eval('case'),';  RMS  Error= 
num2str(E),'%']) 
xlabel('Theta  (degrees)') 
ylabel('|Jtheta|') 

%legend('Jtheta  -  dipsphr2.m','Jback:  no  filtering  -  fem2.m') 
figure 

plot(thdeg,real(Jtheta),'k',thdeg,imag(Jtheta)>'ko',thdeg,real(Jback)> . . . 
'k*',thdeg,imag(Jback),'k+') 

title(['Jexact(real:  imag:  o)  &  Jback(real:  *,  imag:  +):  ',eval('case')]); 
xlabel('Theta  (degrees)') 
ylabel('Jtheta  (Amps/m*m)') 
yn  =  input('Save  data  from  this  run?  (Y/N):  ','s'); 
if  ~isempty(yn) 
if  yn  =  y  |  yn  —  'Y1 

clear  eta  Idl  B  a  zO  dth  th  sth  cth  rd  sthd  cthd  Fth  Etd  Erd  cdth  sdth  ... 

Eti  r  Sn  Pnl  In  n  cn  Ba  Hn  DHn  Ets  Br  Dn  Ett  FIs  FP  Fes  E 
clear  H HI  HH Ha  Hb  He  Hd  JT  Jt  R  R1  R2  alpha  alphal  alphad  an  cp  dphi... 

dtheta  dzf  rb  fho  rhof  rhop  rs  s  talpha  yn  z  zb  zf  zp  zs 
clear  Hpt  Nmax  Nth  f  phi  phipoints... 

radius  rho  theta  thetaFP  thetaFPd  thetad 
save  (eval('case')) 

%save  fO  lambda  al  zOl  rl  Nth  N  inc  ppoints  Hs  Hps  numring  ppoints  numFP... 

%  ppwl  name  thdeg  Hps  Jback  Hphi  Hphitot 

end 

end 
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function  [Fc_s,Fl_s]  =  fcsfl s2(cp,dphi,R) 

%  For  use  with  HPHEERR2.M 
% 

%  This  function  does  the  integration  to  find  Fc(s)  &  Fl(s) 

%  updated  4/19/79  at  home  Dan  Wawrzyniak 

lambda  =  1; 

beta  =  2*pi/lambda; 

phase  =  exp(-i*beta.*R); 

same  =  (l+i*beta.*R)./R.A3; 

FI  =  same.*phase; 

Fl_s  =  (2*dphi*(trapz(Fl 

Fc  =  ones(size(Fl_s))*cp.*same.*phase; 

Fc_s  =  (2  *  dphi  *  (trapz(Fc. '))) . 
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%  DipSpr2.M  computes  scattered  fields  due  to  an  elemental  z-axis 
%  dipole  positioned  above  the  north  pole  of  a  metallic  sphere. 

%  Mod-2  by  M.A.  Morgan  3/14/97 
%  shan3.m  &  legpol2.m  loaded  5/2/97 
clear  all; 

eta=120*pi;  %  free-space  Z_0 

Idl=  1 ;  %  dipole  moment 

fD=300;  %  300  MHz 

B=2*pi*f0/300;  %  wavenumber 

zO=inputCEnter  z-position  of  elemental  dipole  in  meters: '); 
a=inputCEnter  metal  sphere  radius  in  meters: '); 

Nth=input(Enter  number  of  theta  segments: '); 
dth=pi/Nth;  th=(0:dth:pi)';  thdeg=180*th/pi; 
sth=sin(th);  cth=cos(th); 

%  computing  rd,  sin(thd),  and  cos(thd)  for  sphere  surface  (r=a) 
rd=sqrt(a*a+z0*z0-2*a*z0*cth);  sthd=a*sth./rd;  cthd=(a* cth-zO) . /rd; 
%  computing  field  components  in  dipole  centered  coordinates 
Fth=Idl*exp(-j*B*rd)./(4*pi*rd); 

Hpi=Fth.  *  sthd.  *  (j  *B  +  l./rd); 

Etd=eta*Fth.*sthd.*(j*B  +  l./rd  -  j./(B*rd.*rd)); 

Erd=2*eta*Fth.  *  cthd.  *(  1  ./rd  -  j./(B*rd.  *rd)); 

%  translating  incident  field  spherical  components 
cdth=cth.*cthd  +  sth.*sthd;  sdth=sth.*cthd  -  cth.*sthd; 

Eti=Etd.*cdth  -  Erd.*sdth;  %  note  that  Hpi=Hpd 
%  Computing  spherical  harmonic  coefficients  for  scattered  field 
%  Ets  expansion  which  minimizes  the  total  sphere  surface 
%  tangential  field  =  Eti+Ets 
N=input(Enter  number  of  spherical  harmonics: '); 

Sn=-dth*Eti.*sth;  Pnl=legpol2(cth,N);  In=(Pnl.')*Sn; 
n=(l  :N)';  cn=(2*n+l)./(2*n.*(n+l)); 

Ba=B*a;  [Hn,DHn]=shan3(Ba,N); 

Hn  =  Hn.';  DHn  =  DHn.'; 

Dn=DHn/Ba;  an=(cn.*In)./(j*eta*Dn); 

Ets=j*eta*Pnl*(an.*DHn)/Ba;  Ett=Eti+Ets;  %  ideally  Ett  -->  0 

Hps=Pnl*(an.*Hn)/Ba;  Hpt=Hpi+Hps;  %  Jtheta=Hjpt 

Jtheta  =  -Hpt; 

save  dipsphr2  Jtheta  thdeg 

%  Plotting  expansions  on  sphere  surface 

elf  reset;  plot(thdeg,abs(Eti), thdeg, abs(Ett),'.'); 

xlabel('theta  (deg)'); 

ylabel('|E_theta|:  Incident  (solid);  Total  (dots)'); 

title(['Dipole  at  z0=',num2str(z0),';  Sphere  Radius=',num2str(a),... 

';  No.  Modes-, int2str(N)]); 
figure(l); 

hcpy=input('Enter  1  for  Hard  Copy: '); 
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if  hcpy  —  1,  print;  end; 

elf  reset;  plot(thdeg,abs(Hpi),thdeg,abs(Hps),'.',thdeg,abs(Hpt), '--'); 
xlabel('theta  (deg)'); 

ylabel('|H_phi|:  Incident  (solid);  Scattered  (dots);  Total  (dashed)'); 
title([T)ipole  at  zO-,num2str(zO),';  Sphere  Radius- ,num2str(a),... 

';  No.  Modes=',int2str(N)]); 
figure(l); 

hcpy=input(Enter  1  for  Hard  Copy: '); 
if  hcpy  =  1,  print;  end; 

%  Now  compute  scattered  fields  at  specified  radius 
r=inputCEnter  radius  in  meters  to  field  point: '); 

Br=B*r;  [Hn,DHn]=shan3(Br,N); 

Hn  =  Hn.';  DHn  =  DHn.'; 

Ets=j*eta*Pnl  *(an.*DHn)/Br;  Hps=Pnl*(an.*Hn)/Br; 

%  Plotting  scattered  field  expansions 

elf  reset;  plot(thdeg,abs(Ets),thdeg,eta*abs(Hps),'.'); 

xlabel('theta  (deg)'); 

ylabel('|E_theta|  (solid);  eta*|H_phi|  (dots)'); 
title(['Scattered  Fields  at  r=',num2str(r),'  m;  Dipole  at  z0=',... 
num2str(z0),';  Sphere  Radius- ,num2str(a),';  No.  Modes— ,int2str(N)]); 
figure(l); 

hcpy=input('Enter  1  for  Hard  Copy: '); 

if  hcpy  =  1,  print;  end; 

end; 
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%  FEM2.M  provides  finite-element  solution  for  user-specified 
%  metallic  body  of  revolution  using  VARELA  algorithm. 

%  By  M.A.  Morgan  3/20/97-  4/17/97 
clear  all 

%  Generating  surface  and  boundary  points  then  mesh  topology 
sgen=input(Enter  name  of  surface  generation  program:  ','s'); 
eval(['[rs  zs  rb  zb]- ,sgen,';']); 

[rho  zee  nd  elnd]=mesh2(rs,zs,rb,zb);  %  Computing  mesh 
%  Computing  Mesh  Parameters 
N=length(rho);  %  Number  of  Nodes  in  Mesh 

L=length(elnd(:,  1 ));  %  Number  of  Elements  in  Mesh 

Ns=length(nd);  %  No.  of  Surface  or  Boundary  Nodes 

Nm=N-nd(l)-nd(Ns)-Ns+2;  %  Number  of  Nodes  for  Unknown  Hphi 
disp(* ') 

disp(Tinite  Element  Mesh  Parameters: ') 
disp(['  Number  of  Nodes:  ',int2str(N)]) 
disp(['  Number  of  Elements:  ',int2str(L)]) 
disp(['  Number  of  Unknowns:  ',int2str(Nm)]) 
disp(' ') 

dispCPress  a  Key  to  Display  Mesh  or  <Ctl>  C  to  Abort'); 
pause 

%  Plotting  Mesh  Nodes  and  Elements 
%clf  reset 

figure('PaperPosition', [1.5  0.5  5  3.75]) 
if  Ns  <=19,  plot(rho,zee,'go'); 

else,  plot(rho,zee,'k.');  end 
hold  on 
for  1=1  :L 
nds=elnd(l,:); 

rl=[rho(nds);  rho(nds(l))];  zl=[zee(nds);  zee(nds(l))]; 
plot(rl,zl,'k'); 
end 

v=axis;  v(l)=0;  v(2)=v(4)-v(3);  axis(v);  axis  square; 
title('2-D  Axisymmetric  Shape  with  mesh') 
xlabel('meters') 
ylabel('meters') 

text(2.7*v(2)/5,5*v(4)/6,'Mesh  Parameters: ') 

text(3  * v(2)/5 ,3  *  v(4)/4,  [Nodes :  ',int2str(N)]) 

text(3  * v(2)/ 5,2*  v(4)/3 ,  ['Elements :  ',int2str(L)]) 

text(3*v(2)/5,7*v(4)/12, ['Unknowns:  ',int2str(Nm)]) 

yn=inputCPrint  Hard  Copy  ?  (Y/N):  ','s'); 

if  ~isempty(yn),  if  yn  ==  'y'  |  yn  =  'Y',  print;  end;  end 

%  Loading  System  Matrices  By  Indexing  Through  Elements 

disp('Sparse  Matrix  Loading  Using  Element  Contributions  Can  Now  Begin') 

f=input(Enter  Frequency  in  MHz  to  Start  Loading;  <Ctl>  C  to  Stop: '); 
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k2=(pi*f/150)A2;  T=zeros(3,3); 

A=sparse(Nm,Nm);  B=sparse(Nm,Ns-2); 

12=0;  nf=cumsum(nd);  ns=nf-nd+l ;  %  node  finish/start  per  row 

for  i=2:Ns; 

11=12+1;  12=ll+nd(i)+nd(i-l)-3;  %  start/finish  elements 
for  1=11:12 

nds=elnd(l,:);  rp=rho(nds);  zp=zee(nds); 

T=varela(k2,rp,zp); 

dispflToading  Element  ',int2str(l),'  of ',int2str(L)]) 

%  Loading  sparse  matrix  contributions 
for  m=l:3; 
ma=0; 

if  nds(m)  nf(i-l)  &  nds(m)  ~=nf(i),  %  m  not  on  boundary 
if  nds(m)  >  nf(i-l)  &  i  <  Ns, 
ma=nds(m)-ns(2)-i+3;  end 
if  nds(m)  <  ns(i)  &  i  >  2, 
ma=nds(m)-ns(2)-i+4;  end 
if  ma  >  0,  %  m  is  a  solution  node 

for  n=l:3; 
na=0; 

if  nds(n)  ~=  nf(i-l)  &  nds(n)  ~=nf(i), 
if  nds(n)  >  nf(i-l)  &  i  <  Ns, 
na=nds(n)-ns(2)-i+3;  end 
if  nds(n)  <  ns(i)  &  i  >  2, 
na=nds(n)-ns(2)-i+4;  end 
if  na  >  0, 

A(ma,na)=  A(ma,na)+T(m,n);  end 
end  %  n  not  on  boundary  end 
if  nds(n)  =  nf(i-l)  &  i  >  2, 

B(ma,i-2)=B(ma,i-2)-T(m,n);  end 
if  nds(n)  =  nf(i)  &  i  <  Ns, 

B(ma,i-l)=B(ma,i-l)-T(m,n);  end 
end  %  n-loop  end 
end  %  ma  >  0  end 
end  %  m  not  on  boundary  end 
end  %  m-loop  end 
end  %  1-loop  end 
end  %  i-loop  end 

dispCMatrices  Loaded,  Now  Solving  System  ...  Be  Patient') 

%  Solve  sparse  matrix  system  to  construct  transfer  function 
%  array  relating  internal  Hphi  to  Ns-2  nodal  boundary  values 
H  =  A\B;  %  clear  A  B  C  D  %  freeing  memory 
%  Extract  rows  to  form  transfer  array  relating  Ns  nodal  values 
%  of  Hphi  (including  two  z-axis  Hphi=0  BC's)  to  Ns  nodal  values 
%  of  Hphi  on  the  PEC  surface.  Note  that  the  PEC  surface  output 
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%  column  vector  has  Hphi(l)=Hphi(Ns)=0  since  these  correspond 
%  to  nodes  on  the  z-axis. 

Hs=zeros(Ns,Ns); 

row=l;  Hs(2,2:Ns-l)  =  H(row,:); 

for  n=3:Ns-l; 

row=row+nd(n-l)-l;  Hs(n,2:Ns-l)  =  H(row,:); 
end 

n=l:Ns;  elf;  plot(n,diag(Hs),'r') 

v=axis;  axis([l  Ns  v(3)  v(4)]) 

xlabel('row(n)');  ylabelCReal  diag[Hs]') 

title(['Diag[Hs]  for  ',sgen,';  Ns- ,int2str(Ns),';  N-,int2str(N),... 

Nm- ,int2str(Nm),';  L=',int2str(L)]);  figure(l) 
yn=inputCPrint  Hard  Copy  ?  (Y/N): 
if  -isempty(yn),  if  yn  =  'y'  |  yn  —  'Y,  print;  end;  end 
%  Saving  Needed  Parameters  and  Hs  Transfer  Function  Array 
save  femdat  f  rs  zs  rb  zb  Hs 
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%  Program  FEMCHK4.M  compares  exact  and  FEM2.M  computed 
%  nodal  values  for  spherical  harmonic  single  mode  solutions 
%  of  an  offset  metal  sphere.  By  M.A.  Morgan  4/21/97 
%  Mod-2a  (4/30/97)  uses  MatLab  Bessel  and  Hankel  Functions  in  shan3.m 
%  Mod-3  (5/2/97)  uses  Wronskian  to  simplify  r=a  "exact"  solution 
%  and  computes  and  plots  errors  for  applied  modes 
%  Mod-4  (5/3/97)  uses  real  mode  BC's  PnAl(cos(theta))  on  r=b 
%  which  produces  a  real  internal  solution 
clear  all; 

case  =  inputCWhat  case  are  you  running?  ','s'); 
load  femdat;  %  loading  f  rs  zs  rb  zb  Hs 

k=pi*f/150;  Ns=length(rs); 

d=(zs(l)+zs(Ns))/2  -  (zb(l)+zb(Ns))/2;  %  computing  offset 

a=zs(l)-d;  b=zb(l);  Ra=k*a;  Rb=k*b;  %  reset  sphere  centers 
zbd=zb-d;  r_b=sqrt(rb.A2+zbd.A2);  b=zb(l); 
th=linspace(0,pi,Ns)';  c_a=cos(th);  th_a=linspace(0,180,Ns)'; 
c_b=zbd./r_b;  th_b=l  80*acos(c_b)/pi; 

N=inputCEnter  Upper  Spherical  Harmonic  Mode  to  Check: '); 
[hna,Dhna]=shan3  (Ra,N);  [hnb,Dhnb]=shan3  (Rb,N); 

DJna=real(Dhna);  Jnb=real(hnb); 

DNna=-imag(Dhna);  Nnb=-imag(hnb); 

Dn=(Jnb.*DNna-DJna.*Nnb);  An=ones(Ns,l)*DNna;  Bn=ones(Ns,l)*DJna; 
[Hn,DHn]=shan3(k*r_b,N);  Jn=real(Hn);  Nn=-imag(Hn); 
pna=legpol2(c_a,N);  pnb=legpol2(c_b,N);  pct=zeros(N,l); 
Ha=b*pna./(a*ones(Ns,  l)*Dn);  %  "exact"  r=a 

Hb=b,|‘pnb.*(An.*Jn-Bn.*Nn)./(r_b*Dn);  %  "exact"  r=b 
Hc=Hs*Hb;  %  Computed  solution  on  r=a 

for  n=l  :N; 

pct(n)=100*sqrt((Ha(:,n)-Hc(:,n))'*(Ha(:,n)-Hc(:,n)))/... 

sqrt(Ha(:,n)'*Ha(:,n)); 

end 

[tsegl,tseg2]  =  size(Hs); 
figure('PaperPosition',[1.5  0.5  5  3.75]) 
bar(l:N,pct,'k'),  v=axis;  v(l)=0;  v(2)=N+l;  v(3)=0;  axis(v); 
xlabel(Mode  Index  "n"');  ylabel('Percent'); 

title(['Single  Mode  RMS  Error  for  Centered  PEC  Sphere:  ',eval('case')]) 
yn=input(Print  Hard  Copy  ?  (Y/N):  ','s'); 
if  ~isempty(yn),  if  yn  =  'y'  |  yn  =  'Y',  print;  end;  end 
while  1, 

n=input(['Enter  n  <=  ',int2str(N),... 

'  for  Spherical  Harmonic  (0  to  stop): ']); 
ifisempty(n),  break; 
elseif  n  =  0,  break; 

elseif  n  >  N,  disp('n  exceeds  N');  break;  end 
figure('PaperPosition',[1.5  0.5  5  3.75]) 


90 


plot(th_b,Hb(:,n),'k') 

xlabel('Theta  (deg)');  ylabel('Hphi  (A/m)'); 
title(['Single  Mode  Hphi  at  Outer  Boundary:  ',eval('case'), 
n=',int2str(n)]) 

v=axis;  axis([0  180  v(3)  v(4)]); 

yn=inputCPrint  Hard  Copy  ?  (Y/N): 

if  -isempty(yn),  if  yn  —  'y'  |  yn  ==  'Y,  print;  end;  end 

if  length(th_a)  <38, 

figure('PaperPosition',[1.5  0.5  5  3.75]) 

plot(th_a,Ha(:,n),'k',th_a,Hc(:,n),'k.') 

title([Exact  (line),  FEM  (dots):  ',eval(’case'),... 

';  n- ,int2str(n),';  RMS  Error=',num2str(pct(n)),'%']) 

else, 

figureCPaperPosition',[1.5  0.5  5  3.75]) 
plot(th_a,Ha(:  ,n),'k',th_a,Hc( :  ,n),'k. '); 
title([Exact  (line),  FEM  (dots):  ',eval('case'),... 

';  n=',int2str(n),';  RMS  Error=',num2str(pct(n)),'%']) 
end 

xlabel('Theta  (deg)’);  ylabel('Hphi  (A/m)’); 
v=axis;  axis([0  180  v(3)  v(4)]); 
yn=inputCPrint  Hard  Copy  ?  (Y/N):  ','s'); 
if  ~isempty(yn),  if  yn  —  'y'  |  yn  —  'Y',  print;  end;  end 
end 
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%  hspherel.m 

%  Program  Hsphere.M  uses  spherical  harmonics  to  compute 
%  the  Hs  array  for  a  centered  PEC  sphere.  Saves  in  same 
%  format  as  FEM2.m  and  can  be  used  in  FEMCHKxx.m  programs 
%  By  M.  A.  Morgan  5/2/97 
%  5/10/97  modified  by  Dan  Wawrzyniak 
%  Nmax  =  2*k*b 
clear  all 

f=input('Enter  frequency  in  MHz: ');  k=pi*f/150;  lambda  =  3e2/f; 
a=input(Enter  metal  sphere  radius  "a"  in  meters: '); 
b=input('Enter  outer  mesh  radius  "b"  in  meters: '); 

Nsl=inputCEnter  number  of  theta  points  per  wavelength: '); 

Ns  =  fix(pi*b*Nsl/lambda); 
disp(['Total  theta  points  =  ',int2str(Ns)]) 

th=linspace(0,pi,Ns);  thd=linspace(0,180,Ns);  dthl2=pi/(12*(Ns-l)); 
Ra=k*a;  Rb=k*b;  cth=cos(th);  sth=sin(th);  Nmax=2*fix(Rb); 
disp(' ');  disp(['Estimated  Nmax=integer[2*k*b]  is  ',int2str(Nmax)]) 
N=inputCEnter  Upper  Mode  Order  to  Use  in  Computing  Hs: '); 
disp(' ') 

[Hna,DHna]=shan3  (Ra,N);  [Hnb,DHnb]=shan3  (Rb,N); 

Pn  1  =legpol2(cth,N);  n=l:N;  Cn=(2*n+l)./(2*n.*(n+l)); 
Qn=(b/a)./(real(DHna).  *imag(Hnb)-real(Hnb).  *imag(DHna)); 
n=l  :N;  A=zeros(N,Ns);  Hs=zeros(Ns,Ns); 

%  Numerical  integration  of  up(th)*Pnl(cth)*sth 
for  p=2:Ns-l; 

fn=Pnl(p-l:p+l,:);  g=sth(p-l:p+l); 
Inp=dthl2*(fh(l,:)*g(l)+6*fh(2,:)*g(2)+fh(3,:)*g(3)+... 

fii(l,:)*g(2)+fh(2,:)*g(l)+fh(2,:)*g(3)+fh(3,:)*g(2)); 

A(:,p)=(Cn.*Qn.*Inp).'; 

end 

Hs=Pnl*A; 

n=l:Ns;  elf;  plot(n,diag(Hs),'r') 
v=axis;  axis([l  Ns  v(3)  v(4)]) 
xlabel('row(n)');  ylabel('Real  Hs  Diagonal') 
title(['Diag[Hs]  using  Hsphere.m  for  Ns- ,int2str(Ns),... 

Nmax- ,int2str(N),]); 
figure(l) 

yn=inputCPrint  Hard  Copy  ?  (Y/N): 
if  ~isempty(yn),  if  yn  —  'y'  |  yn  =  'Y',  print;  end;  end 
rs=a*sth';  zs=a*cth';  rb=b*sth';  zb=b*cth'; 

%  Saving  Needed  Parameters  and  Hs  Transfer  Function  Array 
save  femdat  f  rs  zs  rb  zb  Hs 
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function  [Hn,DHn]  =  shan3(X,N) 

%  Riccati  Spherical  Hankel  Function  HnA(2)(X)  =  Jn(X)  -j*Yn(X) 

%  and  Derivative.  Using  MatLab  cylindrical  functions  of  order  n+1/2. 
%  Function  returns:  Hn(n)=HnA(2)(x)  and  DHn(n)=dHnA(2)(x)/ dx 
%  for  n=l  to  N.  Input  N  is  a  positive  integer.  X  is  an  array 
%  of  real  or  complex  values  [xl  x2  ...  xM],  Output  Hn(xm)  and 
%  DHn(xm)  are  M  by  N  complex  arrays.  By  M.  A.  Morgan. 

%  Mod-3  (5/1/97)  allows  X  to  be  either  row  or  column  array. 

[ml  m2]=size(X);  M=max(ml,m2);  Hn=zeros(M,N);  DHn=Hn; 
if  ml=l,  x=X';  else,  x=X;  end 
Nu=(l:N)+0.5;  xn=x*ones(l,N); 

J0=sqrt(pi*x/2).  *besselj(0. 5,x); 

Jn=sqrt(pi*xn/2).  *besselj(Nu,x); 

N0=sqrt(pi*x/2).  *bessely(0. 5,x); 

Nn=sqrt(pi*xn/2).*bessely(Nu,x); 

Hn=Jn-j*Nn;  H0=J0-j*N0;  DHn(:,l)=HO-Hn(:,l)./x; 
ifN  >  1,  DHn(:,2)=Hn(:,l)-2*Hn(:,2)./x;  end 
if  N  >  2, 

for  n=3:N;  DHn(:,n)=Hn(:,n-l)-n*Hn(:,n)./x;  end 
end 
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function  Pnl  =  legpol2(X,N) 

%  Computing  a  matrix  of  associated  Legendre  polynomials 
%  PnAm(x)  with  m=l,  for  n=l,  2, ..  N  where  N  is  a  positive  integer 
%  and  X  is  an  M-element  array  of  real  values  [xl  x2  ...  xM], 

%  Computed  real  array  Pnl(xm)  has  M  rows  and  N  columns. 

%  Using  upward  recurrence  formula  from  page  953  of  Balanis  - 

%  Advanced  EM  Engineering,  Wiley,  1989.  Program  by  M.  A.  Morgan 

%  Mod-2  (5/1/97)  allows  either  row  or  column  input  X  array 

[ml  m2]=size(X);  M=max(ml,m2);  Pnl=zeros(M,N); 

if  m  1=1,  x=X';  else,  x=X;  end 

Pnl=zeros(M,N);  Pnl  (:,  l)=-sqrt(l-x.  *x); 

ifN  >  1,  Pnl(:,2)=3*x.*Pnl(:,l);  end 

if  N  >  2, 

forn=2:N-l; 

nl=l/n;  Pnl(:,n+l)=(2+nl)*x.*Pnl(:,nHl+nl)*Pnl(:,n-l); 

end 

end 
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function  [rs,zs,rb,zb]  =  osphere 

%  Computing  meridian  surface  coordinates  for  offset  sphere 
%  (rs,zs)  with  a=radius  and  d=z-axis  offset.  Outer 
%  mesh  boundary  is  centered  b=radius  sphere  (rb,zb). 

%  Returning  coordinates  in  Ns  x  4  array. 

%  By  M.  A.  Morgan  3/24/97 

disp('Offset  Sphere  Surface  and  Mesh  Boundary  Program'); 
a=input('Enter  sphere  radius  (meters): '); 
d=input('Enter  sphere  offset  (meters): '); 
b=input('Enter  mesh  radius  (meters): '); 

Ns=input(Enter  number  of  surface  points: '); 
tha=linspace(0,pi,Ns)';  sta=sin(tha);  cta=cos(tha); 
rs=a*sta;  zs=a*cta+d; 

thb=tha-asin(d*sin(pi-tha)/b);  stb=sin(thb);  ctb=cos(thb); 
rb=b*stb;  zb=b*ctb; 

elf  reset;  plot(rs,zs,rs,zs,'.g',rb,zb,rb,zb,'.g') 
v=axis;  v(l)=0;  v(2)=v(4)-v(3);  axis(v);  axis  square; 
xlabelCPress  a  Key  to  Display  Mesh  or  <Ctrl>  C  to  Abort') 
figure(l) 


function  [rs,zs,rb,zb]  =  cone 

%  Computing  meridian  surface  coordinates  for  centered  cone 
%  (rs,zs)  with  h=height  and  b=base  radius .  Outer 
%  mesh  boundary  is  centered  a=radius  sphere  (rb,zb). 

%  Returning  coordinates  in  Ns  x  4  array. 

%  By  M.  A.  Morgan  3/27/97 

disp('Conical  Surface  and  Spherical  Mesh  Boundary  Program'); 
h=inputCEnter  cone  height  (meters): '); 
b=inputOEnter  cone  base  radius  (meters): '); 
a=inputCEnter  outer  mesh  spherical  radius  (meters): '); 
Ns=inputCEnter  number  of  surface  points: '); 
rs=zeros(Ns,l);  zs=rs;  th=rs; 

S=b+sqrt(b*b+h*h);  %  total  surface  length  on  cone 
Nb=fix(Ns*b/S)+l;  ifNb>Ns-2,  Nb=Ns-2;  end;  db=b/Nb;  Nz=Ns-Nb; 
rs(l  :Nz)=linspace(0,b,Nz);  zs(l  :Nz)=linspace(h/2,-h/2,Nz); 
rs(Nz+l  :Ns)=linspace(b-db,0,Nb);  zs(Nz+l  :Ns)=(-h/2)*ones(Nb,  1); 
%  distort  boundary  node  positions  to  improve  mesh  near  comer 
a2=atan(h/b)/2;  gm=a2+pi/2;  c=(h/2)-b*tan(a2); 
dl=asin(c*  sin(gm)/ a);  thc=gm+dl;  dtb=(pi-thc)/Nb; 
th(l:Nz)=linspace(0,thc,Nz);  th(Nz+l  :Ns)=linspace(thc+dtb,pi,Nb); 
st=sin(th);  ct=cos(th);  rb=a*st;  zb=a*ct; 
elf  reset;  plot^SjZS^zs/r.'.rbjZbjrb.zb/g.') 
v=axis;  v(l)=0;  v(2)=v(4)-v(3);  axis(v);  axis  square;  figure(l) 
xlabel('Press  a  Key  to  Display  Mesh  or  <Ctrl>  C  to  Abort') 
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